#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Heiko Theißen.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##


#############################################################################
##
#M  CanEasilyTestMembership( <permgroup> )
##
InstallTrueMethod(CanEasilyTestMembership,IsPermGroup);


#############################################################################
##
#M  CanEasilyComputePcgs( <permgroup> )
##
##  solvable permgroups can compute a pcgs. (The group may have been told it
##  is solvable without actually doing a solvability test, so we need the
##  implication.)
InstallTrueMethod(CanEasilyComputePcgs,IsPermGroup and IsSolvableGroup);

#M  CanComputeSizeAnySubgroup
InstallTrueMethod(CanComputeSizeAnySubgroup,IsPermGroup);

#############################################################################
##
#M  AsSubgroup( <G>, <U> )  . . . . . . . . . . . .  with stab chain transfer
##
InstallMethod( AsSubgroup,"perm groups",
    IsIdenticalObj, [ IsPermGroup, IsPermGroup ], 0,
    function( G, U )
    local S;
    # test if the parent is already alright
    if HasParent(U) and IsIdenticalObj(Parent(U),G) then
      return U;
    fi;

    if not IsSubset( G, U ) then
      return fail;
    fi;
    S:= SubgroupNC( G, GeneratorsOfGroup( U ) );
    UseIsomorphismRelation( U, S );
    UseSubsetRelation( U, S );
    if HasStabChainMutable( U )  then
        SetStabChainMutable( S, StabChainMutable( U ) );
    fi;
    return S;
end );

#############################################################################
##
#M  CanEasilyComputeWithIndependentGensAbelianGroup( <permgroup> )
##
InstallTrueMethod(CanEasilyComputeWithIndependentGensAbelianGroup,
    IsPermGroup and IsAbelian);

#############################################################################
##
#F  IndependentGeneratorsAbelianPPermGroup( <P>, <p> )  . . . nice generators
##
InstallGlobalFunction( IndependentGeneratorsAbelianPPermGroup,
    function ( P, p )
    local   inds,       # independent generators, result
            pows,       # their powers
            base,       # the base of the vectorspace
            size,       # the size of the group generated by <inds>
            orbs,       # orbits
            trns,       # transversal
            gens,       # remaining generators
            one,        # identity of `P'
            gens2,      # remaining generators for next round
            exp,        # exponent of <P>
            g,          # one generator from <gens>
            h,          # its power
            b,          # basepoint
            c,          # other point in orbit
            i, j, k;    # loop variables

    # initialize the list of independent generators
    inds := [];
    pows := [];
    base := [];
    size := 1;
    orbs := [];
    trns := [];

    # gens are the generators for the remaining group
    gens := GeneratorsOfGroup( P );
    one:= One( P );

    # loop over the exponents
    exp := Maximum( List( gens, g -> LogInt( Order( g ), p ) ) );
    for i  in [exp,exp-1..1]  do

        # loop over the remaining generators
        gens2 := [];
        for j  in [1..Length(gens)]  do
            g := gens[j];
            h := g ^ (p^(i-1));

            # reduce <g> and <h>
            while h <> h^0
              and IsBound(trns[SmallestMovedPoint(h)^h])
            do
                g := g / pows[ trns[SmallestMovedPoint(h)^h] ];
                h := h / base[ trns[SmallestMovedPoint(h)^h] ];
            od;

            # if this is linear independent, add it to the generators
            if h <> h^0  then
                Add( inds, g );
                Add( pows, g );
                Add( base, h );
                size := size * p^i;
                b := SmallestMovedPoint(h);
                if not IsBound( orbs[b] )  then
                    orbs[b] := [ b ];
                    trns[b] := [ one ];
                fi;
                for c  in ShallowCopy(orbs[b])  do
                    for k  in [1..p-1]  do
                        Add( orbs[b], c ^ (h^k) );
                        trns[c^(h^k)] := Length(base);
                    od;
                od;

            # otherwise reduce and add to gens2
            else
                Add( gens2, g );

            fi;

        od;

        # prepare for the next round
        gens := gens2;
        pows := List( pows,i->i^ p );

    od;

    # return the independent generators
    return inds;
end );

#############################################################################
##
#M  IndependentGeneratorsOfAbelianGroup( <G> )  . . . . . . . nice generators
##
InstallMethod( IndependentGeneratorsOfAbelianGroup, "for perm group",
        [ IsPermGroup and IsAbelian ],
    function ( G )
    local   inds,       # independent generators, result
            one,        # identity of `G'
            p,          # prime factor of group size
            gens,       # generators of <p>-Sylowsubgroup
            g,          # one generator
            o;          # its order

    # loop over all primes
    inds := [];
    one:= One( G );
    for p  in Union( List( GeneratorsOfGroup( G ),
            g -> Factors(Order(g)) ) )  do
      if p <> 1  then

        # compute the generators for the Sylow <p> subgroup
        gens := [];
        for g  in GeneratorsOfGroup( G )  do
            o := Order(g);
            while o mod p = 0  do o := o / p; od;
            if g^o <> g^0  then Add( gens, g^o );  fi;
        od;

        # append the independent generators for the Sylow <p> subgroup
        Append( inds,
                IndependentGeneratorsAbelianPPermGroup(
                    GroupByGenerators( gens, one ), p ) );

      fi;
    od;

    # Sort the independent generators by increasing order
    SortBy( inds, Order );

    # return the independent generators
    return inds;
end );

#############################################################################
##
#F  OrbitPerms( <gens>, <d> ) . . . . . . . . . . . . . orbit of permutations
##
InstallGlobalFunction( OrbitPerms, function( gens, d )
    local   max,  orb,  new,  pnt,  img,  gen;

        # get the largest point <max> moved by the group <G>
        max := LargestMovedPoint( gens );

        # handle fixpoints
        if not d in [1..max]  then
            return [ d ];
        fi;

        # start with the singleton orbit
        orb := [ d ];
        new := BlistList( [1..max], [1..max] );
        new[d] := false;

        # loop over all points found
        for pnt  in orb  do

            # apply all generators <gen>
            for gen  in gens  do
                img := pnt ^ gen;

                # add the image <img> to the orbit if it is new
                if new[img]  then
                    Add( orb, img );
                    new[img] := false;
                fi;

            od;

        od;
        return orb;
end );

#############################################################################
##
#F  OrbitsPerms( <gens>, <D> )  . . . . . . . . . . .  orbits of permutations
##
InstallGlobalFunction( OrbitsPerms, function( gens, D )
    local   max,  dom,  new,  orbs,  fst,  orb,  pnt,  gen,  img;

        # get the largest point <max> moved by the group <G>
        max := LargestMovedPoint( gens );
        dom := BlistList( [1..max], D );
        new := BlistList( [1..max], [1..max] );
        orbs := [];

        # repeat until the domain is exhausted
        fst := Position( dom, true );
        while fst <> fail  do

            # start with the singleton orbit
            orb := [ fst ];
            new[fst] := false;
            dom[fst] := false;

            # loop over all points found
            for pnt  in orb  do

                # apply all generators <gen>
                for gen  in gens  do
                    img := pnt ^ gen;

                    # add the image <img> to the orbit if it is new
                    if new[img]  then
                        Add( orb, img );
                        new[img] := false;
                        dom[img] := false;
                    fi;

                od;

            od;

            # add the orbit to the list of orbits and take next point
            Add( orbs, orb );
            fst := Position( dom, true, fst );

        od;

        # add the remaining points of <D>, they are fixed
        for pnt  in D do
          if pnt>max then
            Add( orbs, [ pnt ] );
          fi;
        od;

        return orbs;
end );


#############################################################################
##
#M  SmallestMovedPoint( <C> ) . . . . . . .  for a collection of permutations
##
BindGlobal( "SmallestMovedPointPerms", function( C )
    local   min,  m,  gen;

    min := infinity;
    for gen  in C  do
        if not IsOne( gen ) then
            m := SmallestMovedPointPerm( gen );
            if m < min  then
                min := m;
            fi;
        fi;
    od;
    return min;
end );

InstallMethod( SmallestMovedPoint,
    "for a collection of permutations",
    true,
    [ IsPermCollection ], 0,
    SmallestMovedPointPerms );

InstallMethod( SmallestMovedPoint,
    "for an empty list",
    true,
    [ IsList and IsEmpty ], 0,
    SmallestMovedPointPerms );


#############################################################################
##
#M  LargestMovedPoint( <C> )  . . . . . . .  for a collection of permutations
##
BindGlobal( "LargestMovedPointPerms", function( C )
    local   max,  m,  gen;

    max := 0;
    for gen  in C  do
        if not IsOne( gen )  then
            m := LargestMovedPointPerm( gen );
            if max < m  then
                max := m;
            fi;
        fi;
    od;
    return max;
end );

InstallMethod( LargestMovedPoint,
    "for a collection of permutations",
    true,
    [ IsPermCollection ], 0,
    LargestMovedPointPerms );

InstallMethod( LargestMovedPoint,
    "for an empty list",
    true,
    [ IsList and IsEmpty ], 0,
    LargestMovedPointPerms );


#############################################################################
##
#F  MovedPoints( <C> )  . . . . . . . . . .  for a collection of permutations
##
InstallGlobalFunction(MovedPointsPerms,function( C )
    local   mov,  gen,  pnt;

    mov := [  ];
    for gen  in C  do
        if gen <> One( gen )  then
            for pnt  in [ SmallestMovedPoint( gen ) ..
                           LargestMovedPoint( gen ) ]  do
                if pnt ^ gen <> pnt  then
                    mov[ pnt ] := pnt;
                fi;
            od;
        fi;
    od;
    return Set( mov );
end);

InstallMethod( MovedPoints, "for a collection of permutations", true,
    [ IsPermCollection ], 0,
    MovedPointsPerms );

InstallMethod( MovedPoints, "for an empty list", true,
    [ IsList and IsEmpty ], 0,
    MovedPointsPerms );

InstallMethod( MovedPoints, "for a permutation", true, [ IsPerm ], 0,
function(p)
  return MovedPointsPerms([p]);
end);


#############################################################################
##
#M  NrMovedPoints( <C> )  . . . . . . . . .  for a collection of permutations
##
BindGlobal( "NrMovedPointsPerms", function( C )
    local mov, sma, pnt, gen;

    mov := 0;
    sma := SmallestMovedPoint( C );
    if sma = infinity  then
        return 0;
    fi;
    for pnt  in [ sma .. LargestMovedPoint( C ) ]  do
      for gen in C do
        if pnt ^ gen <> pnt then
          mov:= mov + 1;
          break;
        fi;
      od;
    od;
    return mov;
end );

InstallMethod( NrMovedPoints,
    "for a collection of permutations",
    [ IsPermCollection ],
    NrMovedPointsPerms );

InstallMethod( NrMovedPoints,
    "for an empty list",
    [ IsList and IsEmpty ],
    NrMovedPointsPerms );


#############################################################################
##
#M  LargestMovedPoint( <G> )  . . . . . . . . . . . . for a permutation group
##
InstallMethod( LargestMovedPoint,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> LargestMovedPoint( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  SmallestMovedPoint( <G> ) . . . . . . . . . . . . for a permutation group
##
InstallMethod( SmallestMovedPoint,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> SmallestMovedPoint( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  MovedPoints( <G> )  . . . . . . . . . . . . . . . for a permutation group
##
InstallMethod( MovedPoints,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> MovedPoints( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  NrMovedPoints( <G> )  . . . . . . . . . . . . . . for a permutation group
##
InstallMethod( NrMovedPoints,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> NrMovedPoints( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  OrbitsMovedPoints( <G> )
##
InstallMethod( OrbitsMovedPoints, "for a permutation group", [IsPermGroup],
function(G)
local o,i;
  o:=Orbits(G,MovedPoints(G));
  o:=List(o,Set);
  List(o,IsRange); # save memory if long orbits
  Sort(o);
  MakeImmutable(o);
  for i in o do
    SetIsSSortedList(i,true);
  od;
  SetIsSSortedList(o,true);
  return o;
end);


#############################################################################
##
#M  \=( <G>, <H> )  . .  . . . . . . . . .  test if two perm groups are equal
##
InstallMethod( \=, "perm groups", IsIdenticalObj, [ IsPermGroup, IsPermGroup ],
function ( G, H )
  # avoid new stabilizer chains if they are not computed yet
  if HasSize(G) and HasSize(H) and Size(G)<>Size(H) then return false;fi;

  if not (HasGeneratorsOfGroup(G) and HasGeneratorsOfGroup(H)) then
    TryNextMethod();
  fi;

  if LargestMovedPoint(G)<>LargestMovedPoint(H)
    or OrbitsMovedPoints(G)<>OrbitsMovedPoints(H) then
    return false;
  fi;

  if IsEqualSet( GeneratorsOfGroup( G ), GeneratorsOfGroup( H ) ) then
    return true;
  fi;

  # try to use one existing stab chain
  if HasStabChainMutable(H) or HasStabChainImmutable(H) then
    return Size(G)=Size(H) and ForAll(GeneratorsOfGroup(G),gen->gen in H);
  fi;
  return Size(G)=Size(H) and ForAll(GeneratorsOfGroup(H),gen->gen in G);
end);



#############################################################################
##
#M  BaseOfGroup( <G> )
##
InstallMethod( BaseOfGroup,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> BaseStabChain( StabChainMutable( G ) ) );


#############################################################################
##
#M  Size( <G> ) . . . . . . . . . . . . . . . . . . size of permutation group
##
InstallMethod( Size,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> SizeStabChain( StabChainMutable( G ) ) );


#############################################################################
##
#M  Enumerator( <G> ) . . . . . . . . . . . . enumerator of permutation group
##
BindGlobal( "ElementNumber_PermGroup",
    function( enum, pos )
    local   elm,  S,  len, n;

    S := enum!.stabChain;
    elm := S.identity;
    n:= pos;
    pos := pos - 1;
    while Length( S.genlabels ) <> 0  do
        len := Length( S.orbit );
        elm := LeftQuotient( InverseRepresentative
                       ( S, S.orbit[ pos mod len + 1 ] ), elm );
        pos := QuoInt( pos, len );
        S := S.stabilizer;
    od;

    if pos <> 0 then
      Error( "<enum>[", n, "] must have an assigned value" );
    fi;
    return elm;
end );

BindGlobal( "NumberElement_PermGroup",
    function( enum, elm )
    local   pos,  val,  S,  img;

    if not IsPerm( elm ) then
      return fail;
    fi;

    pos := 1;
    val := 1;
    S := enum!.stabChain;
    while Length( S.genlabels ) <> 0  do
        img := S.orbit[ 1 ] ^ elm;
        #pos := pos + val * ( Position( S.orbit, img ) - 1 );
        if IsBound(S.orbitpos[img]) then
          pos := pos + val * S.orbitpos[img];
          #val := val * Length( S.orbit );
          val := val * S.ol;
          elm := elm * InverseRepresentative( S, img );
          S := S.stabilizer;
        else
          return fail;
        fi;
    od;
    if elm <> S.identity  then
        return fail;
    fi;
    return pos;
end );

InstallMethod( Enumerator,
    "for a permutation group",
    [ IsPermGroup ],
function(G)
local e,S,i,p;
    # get a good base
    S:=G;
    p:=[];
    while Size(S)>1 do
      e:=Orbits(S,MovedPoints(S));
      i:=List(e,Length);
      i:=Position(i,Maximum(i));
      Add(p,e[i][1]);
      S:=Stabilizer(S,e[i][1]);
    od;
    Info(InfoGroup,1,"choose base ",p);

    S:=StabChainOp(G,rec(base:=p));
    e:=EnumeratorByFunctions( G, rec(
             ElementNumber := ElementNumber_PermGroup,
             NumberElement := NumberElement_PermGroup,
             Length        := enum -> SizeStabChain( enum!.stabChain ),
             PrintObj      := function( G )
                                  Print( "<enumerator of perm group>" );
                              end,

             stabChain     := S ) );
  while Length(S.genlabels)<>0 do
    S.ol:=Length(S.orbit);
    S.orbitpos:=[];
    for i in [1..Maximum(S.orbit)] do
      p:=Position(S.orbit,i);
      if p<>fail then
        S.orbitpos[i]:=p-1;
      fi;
    od;
    S:=S.stabilizer;
  od;
  return e;
end);

#############################################################################
##
#M  Iterator( <G> ) . . . . . . . . . . . . . . iterator of permutation group
##
InstallMethod( Iterator,
    "for a trivial permutation group",
    [ IsPermGroup and IsTrivial ],
function(G)
    return IteratorList([()]);
end);

InstallMethod( Iterator,
    "for a permutation group",
    [ IsPermGroup ],
function(G)
    return IteratorStabChain(StabChain(G));
end);

#############################################################################
##
#M  Random( <G> ) . . . . . . . . . . . . . . . . . . . . . .  random element
##
InstallMethodWithRandomSource( Random,
    "for a random source and a permutation group",
    [ IsRandomSource, IsPermGroup ], 10,
    function( rs, G )
    local   S,  rnd;

    # go down the stabchain and multiply random representatives
    S := StabChainMutable( G );
    rnd := S.identity;
    while Length( S.genlabels ) <> 0  do
        rnd := LeftQuotient( InverseRepresentative( S,
                       Random( rs, S.orbit ) ), rnd );
        S := S.stabilizer;
    od;

    # return the random element
    return rnd;
end );


#############################################################################
##
#M  <g> in <G>  . . . . . . . . . . . . . . . . . . . . . . . membership test
##
InstallMethod( \in,
    "for a permutation, and a permutation group",
    true,
    [ IsPerm, IsPermGroup and HasGeneratorsOfGroup], 0,
    function( g, G )
    if g = One( G )  or  g in GeneratorsOfGroup( G )  then
        return true;
    else
        G := StabChainMutable( G );
        return SiftedPermutation( G, g ) = G.identity;
    fi;
end );


#############################################################################
##
#M  ClosureGroup( <G>, <gens>, <options> )  . . . . . .  closure with options
##
##  The following function implements the method and may be called even with
##  incomplete (temp) stabilizer chains.
BindGlobal("DoClosurePrmGp",function( G, gens, options )
    local   C,          # closure of < <G>, <obj> >, result
            P, inpar,   # parent of the closure
            g,          # an element of gens
            o,          # order
            newgens,    # new generator list
            chain;      # the stabilizer chain created

    options:=ShallowCopy(options); # options will be overwritten
    # if all generators are in <G>, <G> is the closure
    gens := Filtered( gens, gen -> not gen in G );
    if IsEmpty( gens )  then
        return G;
    fi;

    # is the closure normalizing?
    if Length(gens)=1 and
       ForAll(gens,i->ForAll(GeneratorsOfGroup(G),j->j^i in G)) then
      g:=gens[1];
      if Size(G)>1 then
        o:=First(Difference(DivisorsInt(Order(g)),[1]),j->g^j in G);
        options.limit:=Size(G)*o;
      else
        o:=First(Difference(DivisorsInt(Order(g)),[1]),j->IsOne(g^j));
        options.limit:=o;
      fi;
      options.size:=options.limit;
      #Print(options.limit,"<<\n");
    fi;

    # otherwise decide between random and deterministic methods
    P := Parent( G );
    inpar := IsSubset( P, gens );
    while not inpar  and  not IsIdenticalObj( P, Parent( P ) )  do
        P := Parent( P );
        inpar := IsSubset( P, gens );
    od;
    if inpar  then
        CopyOptionsDefaults( P, options );
    elif not IsBound( options.random )  then
        options.random := DefaultStabChainOptions.random;
    fi;

    # perhaps <G> is normal in <C> with solvable factor group

#AH 5-feb-96: Disabled (see gap-dev discussion).
#    if DefaultStabChainOptions.tryPcgs
#       and ForAll( gens, gen -> ForAll( GeneratorsOfGroup( G ),
#                   g -> g ^ gen in G ) )  then
#        if inpar  then
#            C := SubgroupNC( P,
#                         Concatenation( GeneratorsOfGroup( G ), gens ) );
#        else
#            C := GroupByGenerators
#                 ( Concatenation( GeneratorsOfGroup( G ), gens ) );
#        fi;
#        pcgs:=TryPcgsPermGroup( [ C, G ], false, false, false );
#        if IsPcgs( pcgs )  then
#          chain:=pcgs!.stabChain;
#          if inpar  then  C := GroupStabChain( P, chain, true );
#                    else  C := GroupStabChain( chain );           fi;
#          SetStabChainOptions( C, rec( random := options.random ) );
#
#          UseSubsetRelation( C, G );
#          return C;
#        fi;
#    fi;

    # make the base of G compatible with options.base
    chain := StructuralCopy( StabChainMutable( G ) );
    if IsBound( options.base )  then
        ChangeStabChain( chain, options.base,
                IsBound( options.reduced ) and options.reduced );
    fi;

    if LargestMovedPoint( Concatenation( GeneratorsOfGroup( G ),
               gens ) ) <= 100  then
        options := ShallowCopy( options );
        options.base := BaseStabChain( chain );
        for g  in gens  do
            if SiftedPermutation( chain, g ) <> chain.identity  then
                StabChainStrong( chain, [ g ], options );
            fi;
        od;
    else
        o:=ValueOption("knownClosureSize");
        if IsInt(o) then options.limit:=o;fi;
        chain := ClosureRandomPermGroup( chain, gens, options );
    fi;

    newgens:=Concatenation(GeneratorsOfGroup(G),gens);
    if Length(chain.generators)<=Length(newgens) then
      if inpar  then  C := GroupStabChain( P, chain, true );
                else  C := GroupStabChain( chain );           fi;
    else
      if inpar  then  C := SubgroupNC(P,newgens);
                else  C := Group( newgens,One(G) );           fi;
      SetStabChainMutable(C,chain);
    fi;
    SetStabChainOptions( C, rec( random := options.random ) );

    UseSubsetRelation( C, G );
    return C;
end );

InstallOtherMethod( ClosureGroup,"permgroup, elements, options",
    true,
    [ IsPermGroup, IsList and IsPermCollection, IsRecord ], 0,
    DoClosurePrmGp);

InstallOtherMethod( ClosureGroup, "empty list",true,
        [ IsPermGroup, IsList and IsEmpty ], 0,
ReturnFirst);

InstallMethod( ClosureGroup, "permgroup, element",true,
  [ IsPermGroup, IsPerm ], 0,
function( G, g )
  return DoClosurePrmGp( G, [ g ], rec() );
end );

InstallMethod( ClosureGroup, "permgroup, elements",true,
        [ IsPermGroup, IsList and IsPermCollection ], 0,
function( G, gens )
  return DoClosurePrmGp( G, gens, rec() );
end );

InstallOtherMethod( ClosureGroup, "permgroup, element, options",true,
  [ IsPermGroup, IsPerm, IsRecord ], 0,
function( G, g, options )
  return DoClosurePrmGp( G, [ g ], options );
end );

InstallOtherMethod( ClosureGroup, "permgroup, permgroup, options", true,
  [ IsPermGroup, IsPermGroup, IsRecord ], 0,
function( G, H, options )
  return DoClosurePrmGp( G, GeneratorsOfGroup( H ), options );
end );

InstallOtherMethod( ClosureGroup, "empty list and options",true,
        [ IsPermGroup, IsList and IsEmpty, IsRecord ], 0,
ReturnFirst);

#############################################################################
##
#M  NormalClosureOp( <G>, <U> ) . . . . . . . . . . . . . . . . in perm group
##
BindGlobal("DoNormalClosurePermGroup",function ( G, U )
    local   N,          # normal closure of <U> in <G>, result
            chain,      # stabilizer chain for the result
            rchain,     # restored version of <chain>, for `VerifySGS'
            options,    # options record for stabilizer chain construction
            gensG,      # generators of the group <G>
            genG,       # one generator of the group <G>
            gensN,      # generators of the group <N>
            genN,       # one generator of the group <N>
            cnj,        # conjugated of a generator of <U>
            random,  k, # values measuring randomness of <chain>
            param,  missing,  correct,  result,  i, one;

    # get a set of monoid generators of <G>
    gensG := GeneratorsOfGroup( G );
    one:= One( G );

    # make a copy of the group to be closed
    N := SubgroupNC( G, GeneratorsOfGroup(U) );
    UseIsomorphismRelation(U,N);
    UseSubsetRelation(U,N);
    SetStabChainMutable( N, StabChainMutable( U ) );
    options := ShallowCopy( StabChainOptions( U ) );
    if IsBound( options.random )  then  random := options.random;
                                  else  random := 1000;            fi;
    options.temp   := true;

    # make list of conjugates to be added to N
    repeat
        gensN := [  ];
        for i  in [ 1 .. 10 ]  do
            genG := SCRRandomSubproduct( gensG, One( G ) );
            cnj  := SCRRandomSubproduct( Concatenation(
                        GeneratorsOfGroup( N ), gensN ), One( G ) ) ^ genG;
            if not cnj in N  then
                Add( gensN, cnj );
            fi;
        od;
        if not IsEmpty( gensN )  then
           N := ClosureGroup( N, gensN, options );
        fi;
    until IsEmpty( gensN );

    # Guarantee that all conjugates are in the normal  closure: Loop over all
    # generators of N
    gensN := ShallowCopy( GeneratorsOfGroup( N ) );
    for genN  in gensN  do

        # loop over the generators of G
        for genG  in gensG  do

            # make sure that the conjugated element is in the closure
            cnj := genN ^ genG;
            if not cnj in N  then
                N := ClosureGroup( N, [ cnj ], options );
                Add( gensN, cnj );
            fi;

        od;

    od;

    # Verify the stabilizer chain.
    chain := StabChainMutable( N );
    if not IsBound( chain.orbits )  then
        if IsBound( chain.aux )  then
            chain := SCRRestoredRecord( chain );
        fi;
        result := chain.identity;
    elif random = 1000  then
        missing := chain.missing;
        correct := chain.correct;
        rchain  := SCRRestoredRecord( chain );
        result  := VerifySGS( rchain, missing, correct );
        if not IsPerm(result) then
            repeat
                result := SCRStrongGenTest2(chain,[0,0,1,10/chain.diam,0,0]);
            until result <> one;
        fi;
        chain := rchain;
    else
        k := First([0..14],x->(3/5)^x <= 1-random/1000);
        if IsBound(options.knownBase) then
            param := [k,4,0,0,0,0];
        else
            param := [QuoInt(k,2),4,QuoInt(k+1,2),4,50,5];
        fi;
        if options.random <= 200 then
            param[2] := 2;
            param[4] := 2;
        fi;
        result := SCRStrongGenTest( chain, param, chain.orbits,
                          chain.basesize, chain.base,
                          chain.correct, chain.missing );
        if result = chain.identity  and  not chain.correct  then
            result := SCRStrongGenTest2( chain, param );
        fi;
        chain := SCRRestoredRecord( chain );
    fi;
    SetStabChainMutable( N, chain );
    if result <> chain.identity  then
        N := ClosureGroup( N, [ result ] );
    fi;

    # give N the proper randomness
    StabChainOptions(N).random:=random;

    # return the normal closure
    UseSubsetRelation( N, U );
    return N;
end );

InstallMethod( NormalClosureOp, "subgroup of perm group",
  IsIdenticalObj, [ IsPermGroup, IsPermGroup ], 0,
function(G,U)
  # test applicability
  if not IsSubset(G,U) then
    TryNextMethod();
  fi;
  return DoNormalClosurePermGroup(G,U);
end);

#############################################################################
##
#M  ConjugateGroup( <G>, <g> )  . . . . . . . . . . . .  of permutation group
##
InstallMethod( ConjugateGroup, "<P>, <g>", true, [ IsPermGroup, IsPerm ], 0,
function( G, g )
local   H,  S;

    H := GroupByGenerators( OnTuples( GeneratorsOfGroup( G ), g ), One( G ) );
    if HasStabChainMutable(G) then
      S := EmptyStabChain( [  ], One( H ) );
      ConjugateStabChain( StabChainMutable( G ), S, g, g );
      SetStabChainMutable( H, S );
    elif HasSize(G) then
      SetSize(H,Size(G));
    fi;
    UseIsomorphismRelation( G, H );
    return H;
end );

#############################################################################
##
#M  CommutatorSubgroup( <U>, <V> )  . . . . . . . . . . . . . for perm groups
##
InstallMethod( CommutatorSubgroup, "permgroups", IsIdenticalObj,
        [ IsPermGroup, IsPermGroup ], 0,
    function( U, V )
    local   C,       # the commutator subgroup
            CUV,     # closure of U,V
            doneCUV, # boolean; true if CUV is computed
            u,       # random subproduct of U.generators
            v,       # random subproduct of V.generators
            comm,    # commutator of u,v
            list,    # list of commutators
            i;       # loop variable

    # [ <U>, <V> ] = normal closure of < [ <u>, <v> ] >.
    C := TrivialSubgroup( U );
    doneCUV := false;

    # if there are lot of generators, use random subproducts
    if Length( GeneratorsOfGroup( U ) ) *
       Length( GeneratorsOfGroup( V ) ) > 10  then
        repeat
           list := [];
           for i in [1..10] do
               u := SCRRandomSubproduct( GeneratorsOfGroup( U ), One( U ) );
               v := SCRRandomSubproduct( GeneratorsOfGroup( V ), One( V ) );
               comm := Comm( u, v );
               if not comm in C then
                   Add( list, comm ) ;
               fi;
           od;
           if Length(list) > 0 then
              C := ClosureGroup( C, list, rec( random := 0,
                                                  temp := true ) );
              if not doneCUV then
                  CUV := ClosureGroup( U, V );
                  doneCUV := true;
              fi;
              # force closure test
              StabChainOptions(C).random:=DefaultStabChainOptions.random;
              C := DoNormalClosurePermGroup( CUV, C );
           fi;
        until IsEmpty( list );
    fi;

    # do the deterministic method; it will also check correctness
    list := [];
    for u  in GeneratorsOfGroup( U )  do
        for v  in GeneratorsOfGroup( V )  do
           comm := Comm( u, v );
           if not comm in C then
               Add( list, comm );
           fi;
        od;
    od;
    if not IsEmpty( list )  then
        C := ClosureGroup( C, list, rec( random := 0,
                                           temp := true ) );
        if not doneCUV then
           CUV := ClosureGroup( U, V );
           doneCUV := true;
        fi;
        # force closure test
        StabChainOptions(C).random:=DefaultStabChainOptions.random;
        C := DoNormalClosurePermGroup( CUV, C );
    fi;

    return C;
end );

#############################################################################
##
#M  DerivedSubgroup( <G> )  . . . . . . . . . . . . . .  of permutation group
##
InstallMethod( DerivedSubgroup,"permgrps",true, [ IsPermGroup ], 0,
    function( G )
    local   D,          # derived subgroup of <G>, result
            g, h,       # random subproducts of generators
            comm,       # their commutator
            list,       # list of commutators
            i,j;  # loop variables

    # find the subgroup generated by the commutators of the generators
    D := TrivialSubgroup( G );

    # if there are >4 generators, use random subproducts
    if Length( GeneratorsOfGroup( G ) ) > 4 then
        repeat
            list := [];
            for i in [1..10] do
                g := SCRRandomSubproduct( GeneratorsOfGroup( G ), One( G ) );
                h := SCRRandomSubproduct( GeneratorsOfGroup( G ), One( G ) );
                comm := Comm( g, h );
                if not comm in D then
                   Add( list, comm );
                fi;
            od;
            if Length(list) > 0 then
              D := ClosureGroup(D,list,rec( random := 0,
                                               temp := true ) );
              # force closure test
              if IsBound(StabChainOptions(G).random) then
                StabChainOptions(D).random:=StabChainOptions(G).random;
              else
                StabChainOptions(D).random:=DefaultStabChainOptions.random;
              fi;
               D := DoNormalClosurePermGroup( G, D );
            fi;
        until list = [];
    fi;

    # do the deterministic method; it will also check random result
    list := [];
    for i in [ 2 .. Length( GeneratorsOfGroup( G ) ) ]  do
         for j  in [ 1 .. i - 1 ]  do
             comm := Comm( GeneratorsOfGroup( G )[i],
                           GeneratorsOfGroup( G )[j] );
             if not comm in D then
                 Add( list, comm );
             fi;
         od;
    od;
    if not IsEmpty( list )  then
        D := ClosureGroup(D,list,rec( random := 0,
                                        temp := true ) );
        # give D a proper randomness
        if IsBound(StabChainOptions(G).random) then
          StabChainOptions(D).random:=StabChainOptions(G).random;
        else
          StabChainOptions(D).random:=DefaultStabChainOptions.random;
        fi;

        D := DoNormalClosurePermGroup( G, D );
    fi;

    if HasIsNaturalSymmetricGroup(G) and Size(D)<Size(G) then IsNaturalAlternatingGroup(D);fi;

    return D;
end );

#############################################################################
##
#M  IsSimpleGroup( <G> )  . . . . . . . test if a permutation group is simple
##
##  This  is  a most interesting function.   It  tests whether  a permutation
##  group is  simple  by testing whether the group is  perfect and then  only
##  looking at the size of the group and the degree of a primitive operation.
##  Basically  it uses  the O'Nan--Scott theorem, which gives  a pretty clear
##  description of perfect primitive groups.  This algorithm is described  in
##  William M. Kantor,
##  Finding Composition Factors of Permutation Groups of Degree $n\leq 10^6$,
##  J. Symbolic Computation, 12:517--526, 1991.
##
InstallMethod( IsSimpleGroup,"for permgrp", true, [ IsPermGroup ], 0,
    function ( G )
    local   D,          # operation domain of <G>
            hom,        # transitive constituent or blocks homomorphism
            d,          # degree of <G>
            n, m,       # $d = n^m$
            simple,     # list of orders of simple groups
            transperf,  # list of orders of transitive perfect groups
            s, t;       # loop variables

    # if <G> is the trivial group, it is not simple
    if IsTrivial( G )  then
        return false;
    fi;

    # first find a transitive representation for <G>
    D := Orbit( G, SmallestMovedPoint( G ) );
    if not IsEqualSet( MovedPoints( G ), D )  then
        hom := ActionHomomorphism( G, D,"surjective" );
        if Size( G ) <> Size( Image( hom ) )  then
            return false;
        fi;
        G := Image( hom );
    fi;

    # next find a primitive representation for <G>
    D := Blocks( G, MovedPoints( G ) );
    while Length( D ) <> 1  do
        hom := ActionHomomorphism( G, D, OnSets,"surjective" );
        if Size( G ) <> Size( Image( hom ) )  then
            return false;
        fi;
        G := Image( hom );
        D := Blocks( G, MovedPoints( G ) );
    od;

    # compute the degree $d$ and express it as $d = n^m$
    D := MovedPoints( G );
    d := Length( D );
    n := SmallestRootInt( d );
    m := LogInt( d, n );
    if 10^6 < d  then
        Error("cannot decide whether <G> is simple or not");
    fi;

    # if $G = C_p$, it is simple
    if    IsPrimeInt( Size( G ) )  then
        return true;

    # if $G = A_d$, it is simple (unless $d < 5$)
    elif  IsNaturalAlternatingGroup(G)  then
        return 5 <= d;

    # if $G = S_d$, it is not simple (except $S_2$)
    elif  IsNaturalSymmetricGroup(G) then
        return 2 = d;

    # if $G$ is not perfect, it is not simple (unless $G = C_p$, see above)
    elif  Size( DerivedSubgroup( G ) ) < Size( G )  then
        return false;

    # if $\|G\| = d^2$, it is not simple (Kantor's Lemma 4)
    elif  Size( G ) = d ^ 2  then
        return false;

    # if $d$ is a prime, <G> is simple
    elif  IsPrimeInt( d )  then
        return true;

    # if $G = U(4,2)$, it is simple (operation on 27 points)
    elif  d = 27 and Size( G ) = 25920  then
        return true;

    # if $G = PSL(n,q)$, it is simple (operations on prime power points)
    elif  (  (d =      8 and Size(G) = (7^3-7)/2          )  # PSL(2,7)
          or (d =      9 and Size(G) = (8^3-8)            )  # PSL(2,8)
          or (d =     32 and Size(G) = (31^3-31)/2        )  # PSL(2,31)
          or (d =    121 and Size(G) =        237783237120)  # PSL(5,3)
          or (d =    128 and Size(G) = (127^3-127)/2      )  # PSL(2,127)
          or (d =   8192 and Size(G) = (8191^3-8191)/2    )  # PSL(2,8191)
          or (d = 131072 and Size(G) = (131071^3-131071)/2)  # PSL(2,131071)
          or (d = 524288 and Size(G) = (524287^3-524287)/2)) # PSL(2,524287)
      and IsTransitive( Stabilizer( G, D[1] ), Difference( D, [ D[1] ] ) )
    then
        return true;

    # if $d$ is a prime power, <G> is not simple (except the cases above)
    elif  IsPrimePowerInt( d )  then
        return false;

    # if we don't have at least an $A_5$ acting on the top, <G> is simple
    elif  m < 5  then
        return true;

    # otherwise we must check for some special cases
    else

        # orders of simple subgroups of $S_n$ with primitive normalizer
        simple := [ ,,,,,
          [60,360],,,,                  #  5: A(5), A(6)
          [60,360,1814400],,            # 10: A(5), A(6), A(10)
          [660,7920,95040,239500800],,  # 12: PSL(2,11), M(11), M(12), A(12)
          [1092,43589145600],           # 14: PSL(2,13), A(14)
          [360,2520,20160,653837184000] # 15: A(6), A(7), A(8), A(15)
        ];

        # orders of transitive perfect subgroups of $S_m$
        transperf := [ ,,,,
          [60],                         # 5: A(5)
          [60,360],                     # 6: A(5), A(6)
          [168,2520],                   # 7: PSL(3,2), A(7)
          [168,8168,20160]              # 8: PSL(3,2), AGL(3,2), A(8)
        ];

        # test the special cases (Kantor's Lemma 3)
        for s  in simple[n]  do
            for t  in transperf[m]  do
                if    Size( G ) mod (t * s^m) = 0
                  and (((t * (2*s)^m) mod Size( G ) = 0 and s <> 360)
                    or ((t * (4*s)^m) mod Size( G ) = 0 and s =  360))
                then
                    return false;
                fi;
            od;
        od;

        # otherwise <G> is simple
        return true;

    fi;

end );

#############################################################################
##
#M  IsSolvableGroup( <G> )  . . . . . . . . . . . . . . . .  solvability test
##

# fast test for finding iteratively nontrivial commutators in a permutation group.
# if this goes on too long the group may not be solvable. This can be much faster
# for large degree groups than to use the full pcgs machinery (which uses the same
# idea but builds a pcgs on the way).
BindGlobal("QuickUnsolvabilityTestPerm",function(G)
local som,elvth,fct,gens,new,l,i,j,a,b,bound;
  # a few moved points
  som:=MovedPoints(G);
  if Length(som) = 0 then SetIsTrivial(G,true); return true; fi;
  bound:=Int(LogInt(Length(som)^5,3)/2); #Dixon Bound
  if Length(som)>100 then
    som:=som{List([1..100],x->Random(1, Length(som)))};
  fi;

  elvth:=One(G);
  fct:=function(G)
    elvth:=elvth*Product([1..Random(1,5)],x->Random(GeneratorsOfGroup(G))^Random(-3,3));
    return elvth;
  end;

  gens:=GeneratorsOfGroup(G);
  # find one nontrivial commutator
  new:=fail;
  i:=1;
  while i<=Length(gens) and new=fail do
    j:=i+1;
    while j<=Length(gens) and new=fail do
      if ForAny(som,x->(x^gens[i])^gens[j]<>(x^gens[j])^gens[i]) then
        new:=Comm(gens[i],gens[j]);
      fi;
      j:=j+1;
    od;
    i:=i+1;
  od;
  if new=fail then return fail;fi;

  # now take commutators with conjugates to get down
  i:=1;
  repeat
    gens:=new;
    j:=1;
    new:=fail;
    l:=[gens];
    while j<=10 and new=fail do
      a:=Random(l)^fct(G);
      b:=First(l,b->ForAny(som,x->(x^a)^b<>(x^b)^a));
      if b<>fail then
        new:=Comm(a,b);
      else
        Add(l,a);
      fi;
      j:=j+1;
    od;
    if new=fail then
      return fail;
    fi;
    i:=i+1;
  until i>bound;
  return true;
end);

InstallMethod( IsSolvableGroup,"for permgrp", true, [ IsPermGroup ], 0,
function(G)
local pcgs;
  if QuickUnsolvabilityTestPerm(G)=true then return false;fi;

  pcgs:=TryPcgsPermGroup( G, false, false, true );
  if IsPcgs(pcgs) then
    SetIndicesEANormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
    SetIsPcgsElementaryAbelianSeries(pcgs,true);
    if not HasPcgs(G) then
      SetPcgs(G,pcgs);
    fi;
    if not HasPcgsElementaryAbelianSeries(G) then
      SetPcgsElementaryAbelianSeries(G,pcgs);
    fi;
    return true;
  else
    return false;
  fi;
end);

#############################################################################
##
#M  IsNilpotentGroup( <G> ) . . . . . . . . . . . . . . . . . nilpotency test
##
InstallMethod( IsNilpotentGroup,"for permgrp", true, [ IsPermGroup ], 0,
    G -> IsPcgs(PcgsCentralSeries(G) ) );

#############################################################################
##
#M  PcgsCentralSeries( <G> )
##
InstallOtherMethod( PcgsCentralSeries,"for permgrp", true, [ IsPermGroup ], 0,
function(G)
local pcgs;
  pcgs:=TryPcgsPermGroup( G, true, false, false );
  if IsPcgs(pcgs) then
    SetIndicesCentralNormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
    SetIsPcgsCentralSeries(pcgs,true);
    if not HasPcgs(G) then
      SetPcgs(G,pcgs);
    fi;
  fi;
  return pcgs;
end);

#############################################################################
##
#M  PcgsPCentralSeriesPGroup( <G> )
##
InstallOtherMethod( PcgsPCentralSeriesPGroup,"for permgrp", true, [ IsPermGroup ], 0,
function(G)
local pcgs;
  pcgs:=TryPcgsPermGroup( G, true, true, Factors(Size(G))[1] );
  if IsPcgs(pcgs) then
    SetIndicesPCentralNormalStepsPGroup(pcgs,pcgs!.permpcgsNormalSteps);
    SetIsPcgsPCentralSeriesPGroup(pcgs,true);
  fi;
  if IsPcgs(pcgs) and not HasPcgs(G) then
    SetPcgs(G,pcgs);
  fi;
  return pcgs;
end);

#T AH, 3/26/07: This method triggers a bug with indices. It is not clear
#T form the description why TryPcgs... should also give a drived series. Thus
#T disabled.
# #############################################################################
# ##
# #M  DerivedSeriesOfGroup( <G> ) . . . . . . . . . . . . . . .  derived series
# ##
# InstallMethod( DerivedSeriesOfGroup,"for permgrp", true, [ IsPermGroup ], 0,
#     function( G )
#     local  pcgs,  series;
#
#     if   (not DefaultStabChainOptions.tryPcgs
#        or ( HasIsSolvableGroup( G )  and  not IsSolvableGroup( G ) ))
#        and not (HasIsSolvableGroup(G) and IsSolvableGroup(G)) then
#         TryNextMethod();
#     fi;
#
#     pcgs := TryPcgsPermGroup( G, false, true, false );
#     if not IsPcgs( pcgs )  then
#         TryNextMethod();
#     fi;
#     SetIndicesEANormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
#     series := EANormalSeriesByPcgs( pcgs );
#     if not HasDerivedSubgroup( G )  then
#         if Length( series ) > 1  then  SetDerivedSubgroup( G, series[ 2 ] );
#                                  else  SetDerivedSubgroup( G, G );  fi;
#     fi;
#     return series;
# end );

#############################################################################
##
#M  LowerCentralSeriesOfGroup( <G> )  . . . . . . . . .  lower central series
##
InstallMethod( LowerCentralSeriesOfGroup,"for permgrp", true, [ IsPermGroup ], 0,
    function( G )
    local  pcgs,  series;

    if    not DefaultStabChainOptions.tryPcgs
       or HasIsNilpotentGroup( G )  and  not IsNilpotentGroup( G )
       and not (HasIsNilpotentGroup(G) and IsNilpotentGroup(G)) then
        TryNextMethod();
    fi;

    pcgs := TryPcgsPermGroup( G, true, true, false );
    if not IsPcgs( pcgs )  then
        TryNextMethod();
    fi;
    SetIndicesCentralNormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
    series := CentralNormalSeriesByPcgs( pcgs );
    if not HasDerivedSubgroup( G )  then
        if Length( series ) > 1  then  SetDerivedSubgroup( G, series[ 2 ] );
                                 else  SetDerivedSubgroup( G, G );  fi;
    fi;
    return series;
end );

InstallOtherMethod( ElementaryAbelianSeries, "perm group", true,
 [ IsPermGroup and IsFinite],
  # we want this to come *before* the method for Pcgs computable groups
  RankFilter(IsPermGroup and CanEasilyComputePcgs and IsFinite)
  -RankFilter(IsPermGroup and IsFinite),
function(G)
local pcgs;
  pcgs:=PcgsElementaryAbelianSeries(G);
  if not IsPcgs(pcgs) then
    TryNextMethod();
  fi;
  return EANormalSeriesByPcgs(pcgs);
end);

#InstallOtherMethod( ElementaryAbelianSeries, fam -> IsIdenticalObj
#        ( fam, CollectionsFamily( CollectionsFamily( PermutationsFamily ) ) ),
#        [ IsList ], 0,
#    function( list )
#    local  pcgs;
#
#    if    IsEmpty( list )
#       or not DefaultStabChainOptions.tryPcgs
#       or (      HasIsSolvableGroup( list[ 1 ] )
#            and not IsSolvableGroup( list[ 1 ] ) )  then
#        TryNextMethod();
##    fi;
#
##    pcgs := TryPcgsPermGroup( list, false, false, true );
#    if not IsPcgs( pcgs )  then  return fail;
#                           else  return ElementaryAbelianSeries( pcgs );  fi;
#end );

#############################################################################
##
#M  PCentralSeries( <G>, <p> )  . . . . . . . . . . . . . .  p-central series
##
InstallMethod( PCentralSeriesOp,"for permgrp", true, [ IsPermGroup, IsPosInt ], 0,
function( G, p )
    local  pcgs;

  if Length(Factors(Size(G)))>1 then
    TryNextMethod();
  fi;
  pcgs:=PcgsPCentralSeriesPGroup(G);
  if not IsPcgs(pcgs) then
    TryNextMethod();
  fi;
  return PCentralNormalSeriesByPcgsPGroup(pcgs);
end);

#############################################################################
##
#M  SylowSubgroup( <G>, <p> ) . . . . . . . . . . . . . .  Sylow $p$-subgroup
##
InstallMethod( SylowSubgroupOp,"permutation groups", true,
  [ IsPermGroup, IsPosInt ], 0,
function( G, p )
local   S;
  if not HasIsSolvableGroup(G) and IsSolvableGroup(G) then
    # enforce solvable methods if we detected anew that the group is
    # solvable.
    return SylowSubgroupOp(G,p);
  fi;
  S := SylowSubgroupPermGroup( G, p );
  S := GroupStabChain( G, StabChainMutable( S ) );
  SetIsNilpotentGroup( S, true );
  if Size(S) > 1 then
    SetIsPGroup( S, true );
    SetPrimePGroup( S, p );
    SetHallSubgroup(G, [p], S);
  fi;
  return S;
end );

InstallGlobalFunction( SylowSubgroupPermGroup, function( G, p )
    local   S,          # <p>-Sylow subgroup of <G>, result
            q,          # largest power of <p> dividing the size of <G>
            D,          # domain of operation of <G>
            O,          # one orbit of <G> in this domain
            B,          # blocks of the operation of <G> on <D>
            f,          # action homomorphism of <G> on <O> or <B>
            T,          # <p>-Sylow subgroup in the image of <f>
            g, g2,      # one <p> element of <G>
            C, C2;      # centralizer of <g> in <G>

    # get the size of the <p>-Sylow subgroup
    q := 1;  while Size( G ) mod (q * p) = 0  do q := q * p;  od;

    # handle trivial subgroup
    if   q = 1  then
        return TrivialSubgroup( G );
    fi;
    # go down in stabilizers as long as possible
    S := StabChainMutable( G );
    while Length( S.orbit ) mod p <> 0  do
        S := S.stabilizer;
    od;
    G := GroupStabChain( G, S, true );

    # handle full group
    if q = Size( G )  then
        return G;
    fi;

    # handle <p>-Sylow subgroups of size <p>
    if q = p  then
        repeat g := Random( G );  until Order( g ) mod p = 0;
        g := g ^ (Order( g ) / p);
        return SubgroupNC( G, [ g ] );
    fi;

    # if the group is not transitive work with the transitive constituents
    D := MovedPoints( G );
    if not IsTransitive( G, D )  then
        Info(InfoGroup,1,"PermSylow: transitive");
        S := G;
        D := ShallowCopy( D );
        while q < Size( S )  do
            O := Orbit( S, D[1] );
            f := ActionHomomorphism( S, O,"surjective" );
            T := SylowSubgroupPermGroup( Range( f ), p );
            S := PreImagesSet( f, T );
            SubtractSet( D, O );
        od;
        return S;
    fi;

    # if the group is not primitive work in the image first
    B := Blocks( G, D );
    if Length( B ) <> 1  then
        Info(InfoGroup,1,"PermSylow: blocks");
        f := ActionHomomorphism( G, B, OnSets,"surjective" );
        T := SylowSubgroupPermGroup( Range( f ), p );
        if Size( T ) < Size( Range( f ) )  then
            T := PreImagesSet( f, T );
            S := SylowSubgroupPermGroup( T , p) ;
            return S;
        fi;
    fi;

    # find a <p> element whose centralizer contains a full <p>-Sylow subgroup
    Info(InfoGroup,1,"PermSylow: element");
    repeat g := Random( G );  until Order( g ) mod p = 0;
    g := g ^ (Order( g ) / p);
    C := Centralizer( G, g );
    while GcdInt( q, Size( C ) ) < q  do
        repeat g2 := Random( C );  until Order( g2 ) mod p = 0;
        g2 := g2 ^ (Order( g2 ) / p);
        C2 := Centralizer( G, g2 );
        if GcdInt( q, Size( C ) ) < GcdInt( q, Size( C2 ) )  then
            C := C2;  g := g2;
        fi;
    od;

    # the centralizer acts on the cycles of the <p> element
    B := List( Cycles( g, D ), Set );
    Info(InfoGroup,1,"PermSylow: cycleaction");
    f := ActionHomomorphism( C, B, OnSets,"surjective" );
    T := SylowSubgroupPermGroup( Range( f ), p );
    S := PreImagesSet( f, T );
    return S;

end );

#############################################################################
##
#M  Socle( <G> )  . . . . . . . . . . .  socle of primitive permutation group
##
InstallMethod( Socle,"test primitive", true, [ IsPermGroup ], 0,
    function( G )
    local   Omega,  deg,  coll,  ds,  L,  z,  ord,
            p,  i;

    Omega := MovedPoints( G );

    if not IsPrimitive( G, Omega )  then
        TryNextMethod();
    fi;

    # Affine groups first.
    L := Earns( G, Omega );
    if L <> []  then
        return L[1];
    fi;

    deg := Length( Omega );
    if deg >= 6103515625  then
        TryNextMethod();
    # the normal closure actually seems to be faster than derived series, so
    # this is not really a shortcut
    #elif deg < 12960000  then
    #    shortcut := true;
    #    if deg >= 3125  then
    #        coll := Collected( Factors( deg ) );
    #        d := Gcd( List( coll, c -> c[ 2 ] ) );
    #        if d mod 5 = 0  then
    #            m := 1;
    #            for c  in coll  do
    #                m := m * c[ 1 ] ^ ( c[ 2 ] / d );
    #            od;
    #            if m >= 5  then
    #                shortcut := false;
    #            fi;
    #        fi;
    #    fi;
    #    if shortcut  then
    #        ds := DerivedSeriesOfGroup( G );
    #        return ds[ Length( ds ) ];
    #    fi;
    fi;

    coll := Collected( Factors(Integers, Size( G ) ) );
    if deg < 78125  then
        p := coll[ Length( coll ) ][ 1 ];
    else
        i := Length( coll );
        while coll[ i ][ 2 ] = 1  do
            i := i - 1;
        od;
        p := coll[ i ][ 1 ];
    fi;
    repeat
        repeat
            z := Random( G );
            ord := Order( z );
        until ord mod p = 0;
        z := z ^ ( ord / p );
    until Index( G, Centralizer( G, z ) ) mod p <> 0;

    # immediately add conjugate as this will be needed anyhow and seems to
    # speed up the closure process
    L := NormalClosure( G, SubgroupNC( G, [ z,z^Random(G) ] ) );
    if deg >= 78125  then
        ds := DerivedSeriesOfGroup( L );
        L := ds[ Length( ds ) ];
    fi;
    if IsSemiRegular( L, Omega )  then
        L := ClosureSubgroup( L, Centralizer( G, L ) );
    fi;
    return L;
end );

#############################################################################
##
#M  FrattiniSubgroup( <G> ) . . . . . . . . . . . .  for permutation p-groups
##
InstallMethod( FrattiniSubgroup,"for permgrp", true, [ IsPermGroup ], 0,
    function( G )
    local   p,  l,  k,  i,  j;

    if not IsPGroup( G ) then
        TryNextMethod();
    elif IsTrivial( G ) then
      return G;
    fi;
    p := PrimePGroup( G );
    l := GeneratorsOfGroup( G );
    k := [ l[1]^p ];
    for i  in [ 2 .. Length(l) ]  do
        for j  in [ 1 .. i-1 ]  do
            Add(k, Comm(l[i], l[j]));
        od;
        Add(k, l[i]^p);
    od;
    return SolvableNormalClosurePermGroup( G, k );
end );

#############################################################################
##
#M  ApproximateSuborbitsStabilizerPermGroup(<G>,<pnt>) . . . approximation of
##  the orbits of Stab_G(pnt) on all points of the orbit pnt^G. (As not
##  all schreier generators are used, the results may be the orbits of a
##  subgroup.)
##
InstallGlobalFunction( ApproximateSuborbitsStabilizerPermGroup,
    function(G,punkt)
    local one, orbit, trans, gen, pnt, img, i, changed,
          processStabGen, currPt,currGen, gens,Ggens;

    one:= One( G );
  if HasStabChainMutable( G ) and punkt in StabChainMutable(G).orbit then
# if we already have a stab chain and bother computing the proper
# stabilizer, a trivial orbit algorithm seems best.
    return Concatenation([[punkt]],
              OrbitsDomain(Stabilizer(G,punkt),
                           Difference(MovedPoints(G),[punkt])));
#    orbit:=Difference(StabChainMutable(G).orbit,[punkt]);
#    stab:=Stabilizer(G,punkt));
#    # catch trivial case
#    if Length(stab)=0 then
#      stab:=[ one ];
#    fi;
#    stgp:=1;
#
#    processStabGen:=function()
#      stgp:=stgp+1;
#      if stgp>Length(stab) then
#        StabGenAvailable:=false;
#      fi;
#      return stab[stgp-1];
#    end;

  else
    # compute orbit and transversal
    Ggens:=GeneratorsOfGroup(G);
    orbit := [ punkt ];
    trans := [];
    trans[ punkt ] := one;
    for pnt  in orbit  do
      for gen  in Ggens do
        img:=pnt^gen;
        if not IsBound( trans[ img ] )  then
          Add( orbit, img );
          trans[img ] := gen;
        fi;
      od;
    od;

  fi;

  currPt:=1;
  currGen:=0;

    processStabGen:=function()
    local punkt,img,a,b,gen;
      if Random(1,2)=1 then
        # next one
        currGen:=currGen+1;
        if currGen>Length(Ggens) then
          currGen:=1;
          currPt:=currPt+1;
          if currPt>Length(orbit) then
            return ();
          fi;
        fi;
        a:=currPt;
        b:=currGen;
      else
        a:=Random(1, Length(orbit));
        b:=Random(1, Length(Ggens));
      fi;
      # map a with b
      img:=orbit[a]^Ggens[b];
      if img=orbit[a] then
        return ();
      else
        punkt:=orbit[a];
        gen:=Ggens[b];
        while punkt<>orbit[1] do
          gen:=trans[punkt]*gen;
          punkt:=punkt/trans[punkt];
        od;
        if orbit[1]^gen<>img then Error("EHEH"); fi;
        punkt:=img;
        while punkt<>orbit[1] do
          gen:=gen/trans[punkt];
          punkt:=punkt/trans[punkt];
        od;
        if orbit[1]^gen<>orbit[1] then Error("UHEH"); fi;
      fi;
      return gen;
    end;

  # as the first approximation take a subgroup generated by four nontrivial
  # Schreier generators
  changed:=0;
  gens:=[];
  while changed<=3 and currPt<=Length(orbit) do
    gen:=();
    while IsOne(gen) and currPt<=Length(orbit) do
      gen:=processStabGen();
    od;
    Add(gens,gen);
    changed:=changed+1;
  od;
  if Length(gens)=0 then
    orbit:=List(orbit,i->[i]);
  else
    orbit:=OrbitsPerms(gens,orbit);
  fi;

  # todo: further fuse
  return orbit;
end );

#############################################################################
##
#M  AllBlocks . . . Representatives of all block systems
##
InstallMethod( AllBlocks, "generic", [ IsPermGroup ],
function(g)
local gens, one, dom,DoBlocks,pool,subo;

  DoBlocks:=function(b)
  local bl,bld,i,t,n;
    # find representatives of all potential suborbits
    bld:=List(Filtered(subo,i->not ForAny(b,j->j in i)),i->i[1]);
    bl:=[];
    if not IsPrime(Length(dom)/Length(b)) then
      for i in bld do
        t:=Union(b,[i]);
        n:=Blocks(g,dom,t);
        if Length(n)>1 and #ok durch pool:ForAll(Difference(n[1],t),j->j>i) and
           not n[1] in pool then
          t:=n[1];
          Add(pool,t);
          bl:=Concatenation(bl,[t],DoBlocks(t));
        fi;
      od;
    fi;
    return bl;
  end;

  one:= One( g );
  gens:= Filtered( GeneratorsOfGroup(g), i -> i <> one );
  if Length(gens)=0 then return []; fi;
  subo:=ApproximateSuborbitsStabilizerPermGroup(g,
          Minimum(List(gens,SmallestMovedPoint)));
  dom:=Union(subo);
  if not IsTransitive(g,dom) then
    Error("AllBlocks, must be transitive");
  fi;
  pool:=[];
  return DoBlocks(subo[1]);
end);

#############################################################################
##
#F  SignPermGroup
##
InstallGlobalFunction( SignPermGroup, function(g)
  if ForAll(GeneratorsOfGroup(g),i->SignPerm(i)=1) then
    return 1;
  else
    return -1;
  fi;
end );

BindGlobal( "CreateAllCycleStructures", function(n)
local i,j,l,m;
  l:=[];
  for i in Partitions(n) do
    m:=[];
    for j in i do
      if j>1 then
        if IsBound(m[j-1]) then
          m[j-1]:=m[j-1]+1;
        else
          m[j-1]:=1;
        fi;
      fi;
    od;
    Add(l,m);
  od;
  return l;
end );

#############################################################################
##
#F  CycleStructuresGroup(G)  list of all cyclestructures occ. in a perm group
##
InstallGlobalFunction( CycleStructuresGroup, function(g)
local l,m,i, one;
  l:=CreateAllCycleStructures(Length(MovedPoints(g)));
  m:=List([1..Length(l)-1],i->0);
    one:= One( g );
  for i in ConjugacyClasses(g) do
    if Representative(i) <> one then
      m[Position(l,CycleStructurePerm(Representative(i)))-1]:=1;
    fi;
  od;
  return m;
end );


# return e is a^e=b or false otherwise
InstallGlobalFunction("LogPerm",
function(a,b)
local oa,ob,q,aa,m,e,tst,s,c,p,g;
  oa:=Order(a);
  ob:=Order(b);
  if oa mod ob<>0 then Info(InfoGroup,20,"1"); return false;fi;
  q:=oa/ob;
  aa:=a^q;

  m:=1; # order modulo which power is OK so far
  e:=1;
  while m<ob do
    tst:=aa^e/b;
    s:=SmallestMovedPoint(tst); # a point on which there might be a difference
    if s=infinity then # found it
      return q*e mod oa;
    fi;
    c:=Cycle(aa,s);
    g:=Gcd(m,Length(c));
    p:=Position(c,s^b);
    if p=fail or (e mod g<>(p-1) mod g) then return false;fi;
    e:=ChineseRem([m,Length(c)],[e,p-1]);
    m:=Lcm(m,Length(c));
  od;
  if aa^e<>b then return false;fi;
  return q*e mod oa;
end);


#############################################################################
##
#M  SmallGeneratingSet(<G>)
##
InstallMethod(SmallGeneratingSet,"random and generators subset, randsims",true,
  [IsPermGroup],0,
function (G)
local  i, j, U, gens,a,mintry,orb,orp,isok;

  gens := ShallowCopy(Set(GeneratorsOfGroup(G)));

  # remove obvious redundancies...
  # sort elements by descending order; trivial permutations
  # at the end
  SortBy(gens, g -> -Order(g));

  i := 1;
  while i < Length(gens) do
    j := i + 1;
    while j <= Length(gens) do
      a:=LogPerm(gens[i], gens[j]);
      if a = false then
        j := j + 1;
      else
        Remove(gens, j);
      fi;
    od;
    i := i + 1;
  od;

  if Length(gens)<=Length(AbelianInvariants(G))+2 then
    return gens;
  fi;

  # try pc methods. The solvability test should not exceed cost, nor
  # the number of points.
  if #Length(MovedPoints(G))<50000 and
   #((HasIsSolvableGroup(G) and IsSolvableGroup(G)) or IsAbelian(G))
      IsSolvableGroup(G)
      and Length(gens)>3 then
    return MinimalGeneratingSet(G);
  fi;

  # store orbit data
  orb:=Set(Orbits(G,MovedPoints(G)),Set);
  # longer primitive orbits
  orp:=Filtered([1..Length(orb)],x->Length(orb[x])>100 and IsPrimitive(Action(G,orb[x])));

  isok:=function(U)
    if Length(MovedPoints(G))>Length(MovedPoints(U)) or
      Length(orb)>Length(Orbits(U,MovedPoints(U))) or
        ForAny(orp,x->not IsPrimitive(U,orb[x])) then
      return false;
    fi;

    StabChainOptions(U).random:=100; # randomized size
    return Size(U)=Size(G);
  end;

  mintry:=2;
  if Length(gens)>2 then
    # lower bound on what to accept -- use index of G' instead of full
    # structure of G/G', as the latter is more expensive.
    mintry:=Maximum(List(Collected(Factors(Size(G)/Size(DerivedSubgroup(G)))),x->x[2]));
    mintry:=Maximum(mintry,2);
    if mintry=Length(gens) then return gens;fi;
    i:=Maximum(2,mintry);
    while i<=mintry+1 and i<Length(gens) do
      # try to find a small generating system by random search
      j:=1;
      while j<=5 and i<Length(gens) do
        U:=Subgroup(G,List([1..i],j->Random(G)));
        if isok(U) then
          gens:=Set(GeneratorsOfGroup(U));
        fi;
        j:=j+1;
      od;
      i:=i+1;
    od;
  fi;

  # try to delete generators
  i := 1;
  while i <= Length(gens) and Length(gens)>mintry do
    # random did not improve much, try removing generators one by one
    U:=Subgroup(G,gens{Difference([1..Length(gens)],[i])});
    if isok(U) then
      gens:=Set(GeneratorsOfGroup(U));
    else
      i:=i+1;
    fi;
  od;
  return gens;
end);

#############################################################################
##
#M  GeneratorsSmallest(<G>) . . . . . . . . . . . . . for permutation groups
##
DeclareGlobalFunction( "GeneratorsSmallestStab" );

InstallGlobalFunction( GeneratorsSmallestStab, function ( S )
    local   gens,       # smallest generating system of <S>, result
            gen,        # one generator in <gens>
            orb,        # basic orbit of <S>
            pnt,        # one point in <orb>
            T,          # stabilizer in <S>
            o2,img,oo,og; # private orbit algorithm

    # handle the anchor case
    if Length(S.generators) = 0  then
        return [];
    fi;

    # now get the smallest generating system of the stabilizer
    gens := GeneratorsSmallestStab( S.stabilizer );

    # get the sorted orbit (the basepoint will be the first point)
    orb := Set( S.orbit );
    SubtractSet( orb, [S.orbit[1]] );

    # handle the other points in the orbit
    while Length(orb) <> 0  do

        # take the smallest point (coset) and one representative
        pnt := orb[1];
        gen := S.identity;
        while S.orbit[1] ^ gen <> pnt  do
           gen := LeftQuotient( S.transversal[ pnt / gen ], gen );
        od;

        # the next generator is the smallest element in this coset
        T := S.stabilizer;
        while Length(T.generators) <> 0  do
            pnt := Minimum( OnTuples( T.orbit, gen ) );
            while T.orbit[1] ^ gen <> pnt  do
                gen := LeftQuotient( T.transversal[ pnt / gen ], gen );
            od;
            T := T.stabilizer;
        od;

        # add this generator to the generators list and reduce orbit
        Add( gens, gen );

        #SubtractSet( orb,
        #    Orbit( GroupByGenerators( gens, S.identity ), S.orbit[1] ) );

        # here we want to be really fast, so use a private orbit algorithm
        o2:=[S.orbit[1]];
        RemoveSet(orb,S.orbit[1]);
        for oo in o2 do
          for og in gens do
            img:=oo^og;
            if img in orb then
              Add(o2,img);
              RemoveSet(orb,img);
            fi;
          od;
        od;

    od;

    # return the smallest generating system
    return gens;
end );

InstallMethod(GeneratorsSmallest,"perm group via minimal stab chain",
  [IsPermGroup],
function(G)
  # call the recursive function to do the work
  return GeneratorsSmallestStab(MinimalStabChain(G));
end);

InstallMethod(LargestElementGroup,"perm group via minimal stab chain",
  [IsPermGroup],
  # call the recursive function to do the work
  G -> LargestElementStabChain( MinimalStabChain( G ), One( G ) ) );


#############################################################################
##
#M  KnowsHowToDecompose( <G>, <gens> )
##
InstallMethod(KnowsHowToDecompose,
    "perm group and generators: always true",
    IsIdenticalObj,
    [ IsPermGroup, IsList ],
    ReturnTrue);


#############################################################################
##
#M  ViewObj( <G> )
##
InstallMethod( ViewString, "for a permutation group", true,
    [ IsPermGroup and HasGeneratorsOfGroup ], 0,
function(G)
  local gens,str;
  gens:= GeneratorsOfGroup( G );
  if 30 < Length( gens ) * LargestMovedPoint( G ) / GAPInfo.ViewLength then
    str:="<permutation group";
    if HasSize(G) then
      Append(str," of size ");
      Append(str,String(Size(G)));
    fi;
    Append(str," with ");
    Append(str,Pluralize(Length( gens ), "generator"));
    Append(str,">");
  else
    str:="Group(";
    if IsEmpty( gens ) or ForAll(gens,i->Order(i)=1) then
      Append(str,"()");
    else
      Append(str,ViewString(gens));
    fi;
    Append(str,")");
  fi;
  return str;
end);


#############################################################################
##
#M  AsList( <G> ) elements of perm group
##
InstallMethod( AsList, "permgp: ElementsStabChain", [ IsPermGroup ],
  G -> ElementsStabChain( StabChainMutable(G)) );

#############################################################################
##
#M  AsSSortedList( <G> ) elements of perm group
##
InstallMethod( AsSSortedList, "via stabchain", [ IsPermGroup ],
  AsSSortedListNonstored );

InstallMethod( AsSSortedListNonstored, "via stabchain", [ IsPermGroup ],
function( G )
  local list;
  list := ElementsStabChain( StabChainMutable(G));
  Sort(list);
  return list;
end );

#############################################################################
##
#M  ONanScottType( <G> )
##
InstallMethod(ONanScottType,"primitive permgroups",true,[IsPermGroup],0,
function(G)
local dom,s,cs,t,m,stb;
  dom:=MovedPoints(G);
  if not IsPrimitive(G,dom) then
    Error("<G> must be primitive");
  fi;
  s:=Socle(G);
  if IsAbelian(s) then
    return "1";
  elif IsNonabelianSimpleGroup(s) then
    return "2";
  elif Length(dom)=Size(s) then
    return "5";
  else
    # so the group now is of type 3 or 4. Next we determine a simple socle
    # factor and a direct product of all but one socle factor.
    # simple socle factor.

    # as default stab chain chooes pts lexicographically, this is likely a
    # stored stabilizer
    stb:=Stabilizer(s,SmallestMovedPoint(s));
    cs:=CompositionSeries(s);
    t:=cs[Length(cs)-1];
    m:=cs[2];

    # now test Dixon&Mortimer, p.126, Case 1: R1 (projection of pt. stab
    # onto one direct factor) is proper subgroup of T1.

    #if Size(ClosureSubgroup(m,stb))<Size(s) then
    # Since m is normal, this projection is a proper subgroup iff m does not
    # act transitively
    if not IsTransitive(m,dom) then
      return "4c"; # Product action of type 2 with transitive
    fi;

    #  Now we are in case 2, R1=T1. Group must be diagonal (3) or product
    #  action of diagonal with transitive (4). For diagonal case the point
    #  stabilizer must have the order of t, for product case it must be a
    #  power of the order of t
    if Size(stb)=Size(t) then # type 3
      # in the a case the socle is not minimal
      #if Size(NormalClosure(G,t))<Size(s) then
      if Length(Orbit(G,t))<Length(cs)-1 then
        return "3a";
      else
        return "3b";
      fi;
    else
      # Same argument as for 3:
      #if Size(NormalClosure(G,t))<Size(s) then
      if Length(Orbit(G,t))<Length(cs)-1 then
        return "4a";
      else
        return "4b";
      fi;
    fi;
  fi;
end);

#############################################################################
##
#M  SocleTypePrimitiveGroup( <G> )
##
InstallMethod(SocleTypePrimitiveGroup,"primitive permgroups",true,
  [IsPermGroup],0,function(G)
local s,cs,t,id,r;
  if not IsPrimitive(G) then
    ErrorNoReturn("<G> must be primitive on its moved points");
  fi;
  s:=Socle(G);
  cs:=CompositionSeries(s);
  t:=cs[Length(cs)-1];
  id:=IsomorphismTypeInfoFiniteSimpleGroup(t);
  r:=rec(series:=id.series,
             width:=Length(cs)-1,
             name:=id.name);
  if IsBound(id.parameter) then
    r.parameter:=id.parameter;
  else
    r.parameter:=id.name;
  fi;
  return r;
end);

# this function is needed for the transitive and primitive groups libraries
BindGlobal("STGSelFunc",function(a,b)
  if IsFunction(b) then
    return b(a);
  elif IsList(b) and not IsString(b) then
   if IsList(a) and a=b then
     return true;
   fi;
   return a in b;
 else
   return a=b;
 fi;
end);

InstallGlobalFunction(DiagonalSocleAction,function(g,nr)
local ran,d,emb,u,act;
  ran:=[1..nr];
  d:=DirectProduct(List(ran,i->g));
  emb:=List(ran,i->Embedding(d,i));
  u:=List(GeneratorsOfGroup(g),i->Product(ran,j->Image(emb[j],i)));
  u:=SubgroupNC(d,u);
  act:=ActionHomomorphism(d,RightTransversal(d,u),OnRight,"surjective");
  return Image(act,d);
end);

InstallGlobalFunction(ReducedPermdegree,function(g)
  local dom, orb, gdeg, p, ndom, s,hom;
  dom:=MovedPoints(g);
  orb:=ShallowCopy(Orbits(g,dom));
  if Length(orb)=1 then return fail;fi;
  gdeg:=20*LogInt(Size(g),2); # good degree
  SortBy(orb, Length);
  p:=PositionProperty(orb,i->Length(i)>=gdeg);
  if p=fail then p:=Length(orb);fi;
  ndom:=orb[p];
  if Length(ndom)*2>Length(dom) then return fail;fi;
  s:=Stabilizer(g,orb[p],OnTuples);
  p:=p+1;
  if p>Length(orb) then p:=1;fi;
  while Size(s)>1 do
    while not ForAny(orb[p],j->ForAny(GeneratorsOfGroup(s),k->j^k<>j)) do
      p:=p+1;
      if p>Length(orb) then p:=1;fi;
    od;
    ndom:=Union(ndom,orb[p]);
    s:=Stabilizer(s,orb[p],OnTuples);
  od;
  if Length(ndom)*2>Length(dom) then return fail;fi;
  hom:=ActionHomomorphism(g,ndom,"surjective");
  SetIsInjective(hom,true);
  SetSize(Image(hom),Size(g));
  return hom;
end);
